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The clustering properties of sterile neutrinos are studied within a simple extension of the minimal 

standard model, where these neutrinos are produced via the decay of a gauge singlet scalar. The 

distribution function after decoupling is strongly out of equilibrium and features an enhancement at 

small comoving momentum oc 1/y^. Dark matter abundance and phase space density constraints 

from dwarf spheroidal galaxies constrain the mass in the keV range consistent with a Yukawa 

coupling to a gauge singlet with mass and vacuum expectation value in the range ~ 100 GeV and 

a decoupling temperature of this order. The dark matter transfer function and power spectrum 

. . - are obtained from the solution of the non-relativistic Boltzmann-Vlasov equation in the matter 

OO ' dominated era. The small momentum enhancement of the non-equilibrium distribution function 

'n— ^ , leads to long range memory of gravitational clustering and a substantial enhancement of the power 

>— ^ ■ spectrum at small scales as compared to a thermal relic or sterile neutrino produced via non-resonant 

mixing with active neutrinos. The scale of suppression of the power spectrum for a sterile neutrino 

^ , with m ~ keV produced by scalar decay that decouples at ~ lOOGeV is A ~ 488 kpc. At large 

(D . scales T{k) ~ 1 - C k^ /k},{teq) H with C ~ 0(1). At smaU scales 65kpc < A < 500 kpc 

[] ' corrections to the fluid description and memory of gravitational clustering become important, and 

we find T(k) ~ 1.902 e~*'''°J''^ '*="', where kfs{teq) ~ 0.013/kpc is the free streaming wavevector at 
matter-radiation equality. The enhancement of power at small scales may provide a possible relief 
to the tension between the constraints from X-ray and Lyman-a forest data. 
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I. INTRODUCTION 



In the concordance ACDM standard cosmological model dark matter (DM) is composed of primordial particles which 
are cold and collisionlcss 1] . In this cold dark matter (CDM) scenario structure formation proceeds in a hierarchical 
X^J ] "bottom up" approach: small scales become non-linear and collapse first and their merger and accretion leads to 
v^ ■ structure on larger scales. CDM particles feature negligible small velocity dispersion leading to a power spectrum 
^b ' that favors small scales. In this hierarchical scenario, dense clumps that survive the merger process form satellite 
\^ . galaxies. 
f^ ' Large scale simulations seemingly yield an overprediction of satellite galaxies!^ by almost an order of magnitude 

p [ , larger than the number of satellites that have been observed in Milky- Way sized galaxies 0, y, |j, la, IS]- Simulations 

f~^ ■ within the ACDM paradigm also yield a density profile in virialized (DM) halos that increases monotonically towards 

OO '. the center^, H M, B [ly| and features a cusp, such as the Navarro-Frenk- White (NFW) profileQ or more general 

^D I central density profiles p{r) ~ r~^ with 1 < /3 < 1-50, 0, \l^. These density profiles accurately describe clusters 

^ \ of galaxies but there is an accumulating body of observational evidence [ll|, [l^, Uj, Hj \1M, US, Hj \1m that seem to 

indicate that the central regions of (DM)-dominated dwarf spheroidal satellite (dSphs) galaxies feature smooth cores 

instead of cusps as predicted by (CDM). More recently [1£| a "galaxy size" problem has been reported, where large 

'^ , scale simulations at z = 3 yield galaxies that are too small, this problem has been argued to be related to that of the 

■ - ' missing dwarf galaxies. 

Warm dark matter (WDM) particles were invoked |20l [2ll ^^ as possible solutions to these discrepancies both in 
the over abundance of satellite galaxies and as a mechanism to smooth out the cusped density profiles predicted by 
(CDM) simulations into the cored profiles that fit the observations in (dShps). (WDM) particles feature a range of 
velocity dispersion in between the (CDM) and hot dark matter (HDM) leading to free streaming scales that smooth 
out small scale features and could be consistent with core radii of the (dSphs). If the free streaming scale of these 
particles is smaller than the scale of galaxy clusters, their large scale structure properties are indistinguishable from 
(CDM) but may affect the small scale power spectrumj23| so as to p rovide an explanation of the smoother inner 
profiles of (dSphs), fewer satellites and the size of galaxies at z = 3|19|. 

Sterile neutrinos with masses ~ keV may be suitable (WDM) candidates |24l [2a. \2l\. \2al. [29. [3(j| . The main property 
that is relevant for structure formation of any dark matter candidate is its distribution function after decoupling 33, 
|33|, which depends on the production mechanism and the (quantum) kinetics of its evolution from production to 
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decoupling. There is a variety of mechanisms of sterile neutrino production[2J, [2^, [2y, [23, yj], and mixing between 
sterile and active neutrinos can be one of them[2J, [2^, [2^. There is considerable tension between the X-ray [33| and 
Lyman-a forest [30, [33| data if sterile neutrinos are produced via the Dodelson-Widrow (DW)[24| non-resonant mixing 
mechanism, leading to the suggestion JSq that these cannot be the dominant (DM) component. Constraints from the 
Lyman-a forest spectra are particularly im por tant because of its sensitivity to the suppression of the power spectrum 
by free-streaming in the linear regime [36l |37|. The most recent constraints from the Lyman-a forest [37| improve 
upon previous ones, but rely on the Dodelson-Widrow'24] model for the distribution function of sterile neutrinos, 
leaving open the possibility of evading these tight constraints with non-equilibrium distribution functions from other 
production mechanisms. 

The gravitational clustering properties of coUisionlcss (DM) in the linear regime are described by the power spectrum 
of gravitational perturbations. Free streaming of coUisionlcss (DM) leads to a suppression of the transfer function on 
length scales smaller than the free streaming scale via Landau damping[23, [3^, |40|- This scale is determined by the 
decoupling temperature, the particle's mass and the distribution function at decoupling [4 1|. 

In this article we study the gravitational clustering properties of sterile neutrinos as (WDM) candidates within a 
model in which these are produced via the decay of a gauge singlet scalar. Such model has been advocated recently 
in ref.[2^, [30, [sj, and is based on a phenomenologically appealing extension of the minimal standard model[27| 
consistent with the observed neutrino masses and mixing. This model also shares many features in common with 
models of gravitino production[42], hence it provides a viable extension of the standard model to study in detail 
the production and clustering properties of potential (WDM) candidates. In this model sterile neutrinos decouple a 
temperatures much larger (~ 100 GeV ) than in the Dodelson-Widrow scenario (~ 150MeV)|2J], therefore they are 
colder at matter-radiation equality (and today) , being dubbed, for this reason, "chilled" neutrinos in refs. [2^, [30| ■ 
Clustering properties of active neutrinos in non-standard cosmology, for example quintessence have been reported in 
ref.[3li. 

A program that yields a quantitative assessment of a particle physics candidate for (DM) in the linear regime 
implements the following steps: 

• Establish the quantum kinetic equations that describe the production of these particles and follows the evolution 
of their distribution function through their decoupling from the cosmological plasma. 

• The distribution function after decoupling becomes the unperturbed distribution, which determines the abun- 
dance, the primordial phase space densities [33, [33[ and free streaming lengths j41l]. The generalized Tremaine- 
Gunn[4^ constraints obtained in[33| in combination with the recent photometric observations of the phase space 
densities of (dSphs) combined with the DM abundance lead to bounds on the mass, couplings and decoupling 
temperature 1331 . 

• The unperturbed distribution function is input in the Boltzmann-Vlasov equation for density and gravitational 
perturbations [44]. The solution of which yields the transfer function, and the power spectrum. 

We follow this program in the model of sterile neutrino production proposed in refs.[27|, [29|, [30, [3J]. In principle, in 
order to obtain the transfer function and the power spectrum of density perturbations the coupled set of Boltzmann 
equations for baryons, photons, dark matter and gravitational perturbations must be solved [44 [45|. Photons and 
baryons are coupled by Thompson scattering and dark matter only couples to the gravitational perturbations that are 
sourced by all the components. In practice this is a computationally daunting task because popular codes [46] based 
on the set of coupled Boltzmann equations for photons, baryons and dark matter [44| need to be modified to input 
arbitrary non-equilibrium distribution functions, masses and couplings. 

Recently a simple analytic framework to obtain the dark matter transfer function, and consequently the power spec- 
trum during matter domination has been presented [47|. The main premise of this formulation is that the contribution 
from baryons and photons modifies the DM transfer function at most by a few percent |48l l49l| during matter domina- 
tion and that a preliminary robust assessment of the clustering properties of a DM candidate can be systematically 
established by neglecting in first approximation the contribution from baryons and photons. The influence of baryons 
on the DM power spectrum is more prominent on the scale of baryon acoustic oscillations (BAO), corresponding to 
the scale of the sound horizon at recombination, or ~ 150 Mpc today t5(|. On this scale the DM power spectrum does 
not distinguish between (CDM) or (WDM), and at smaller scales, of interest for the satellite and cusp problems, the 
(BAO) features are not prominent and can be safely neglected. 

The main ingredient to study the (DM) transfer function in absence of baryons is the non-relativistic Boltzmann- 
Vlasov equation for DM density and gravitational perturbations. The non-relativistic limit is warranted for particles 
that decoupled early and became non-relativistic prior to matter-radiation equality and for perturbations that entered 
the horizon prior to matter-radiation equality, these describe all the relevant scales for structure formation. 

The method developed in ref. [47| yields a simple analytic approximation to the transfer function that is remarkably 
accurate in a wide range of scales relevant to structure formation. An important ingredient is a non-local kernel that 



depends on the unperturbed distribution function of the decoupled particles. This kernel describes memory of gravi- 
tational clustering and is a correction to the fluid description which become important at small scales. Distribution 
functions that feature larger support at small momentum yield longer range memory kernels thereby enhancing the 
transfer function at small scales 47] placing greater importance on non-equilibrium aspects of the distribution function. 
In this article we study the clustering properties of sterile neutrinos in the model proposed in references [271 l29l . ISOl. |34| 
by implementing all the steps described above, from obtaining and solving the quantum kinetic equation for production 
that establishes the distribution function at decoupling, narrowing the range of parameters, masses and couplings with 
the observational constraints from (DM) abundance and coarse grained phase space densities of (DM) dominated 
fdSphs)|14l. l3a |. and solving the Boltzmann-Vlasov equation obtaining the transfer function and power spectrum 



which we compare to the case of thermal relics or Dodelson- Widrow-tvpe J24] distribution functions 



Summary of results: 

The production of sterile neutrinos of to ~ keV via the decay of scalar gauge singlet with M ^ 100 GeV leads to 
decoupling at a temperature ^ 100 GeV and a distribution function that is strongly out of equilibrium and behaves 
as l/-v/p for small comoving momenta p. 

The constraints from (DM) abundance and coarse-grained phase space density from the latest compilation of 
photometric data from (dSphs) lead to a narrow window in the keV range for the value of the mass of the sterile 
neutrino, consistently with the phenomenologically motivated extension beyond the standard model studied. 

The (DM) transfer function and power spectrum are obtained from the solution of the non-relativistic Boltzmann- 
Vlasov equation for (DM) density and gravitational perturbations during matter domination. We implement a simple 
analytic approximate method [4 7| to obtain the density and gravitational perturbations and the transfer function that 
is remarkably accurate in a wide range of cosmologically relevant scales, as confirmed by the exact solution. This 
approach yields a wealth of information that relates the small scale behavior of the transfer function to the range of 
memory of gravitational clustering, which is determined by the small (comoving) momentum region of the distribution 
function. 

The enhancement of the non-equilibrium distribution at small momentum leads to a long range memory of grav- 
itational clustering and slower fall off of the free-streaming solution. Both features lead to an enhancement of the 
transfer function and power spectrum at small scales. 

We compare the transfer function and power spectrum from sterile neutrinos produced via gauge singlet decay to 
that of relativistic fermions decoupled in local thermodynamic equilibrium (LTE) (thermal relic) and sterile neutrinos 
produced by non-resonant mixing with active neutrinos a la Dodelson- Widrowf245. Thermal relics and (DW)-produced 
sterile neutrinos feature the same transfer function for similar ratios of the decoupling temperature to mass. 

The transfer function and power spectrum for sterile neutrinos produced by scalar decay is substantially enhanced 
with respect to that of (DW)-sterile neutrinos (and thermal relics) at small scales A < 500 kpc. Whereas for (DW)- 
sterile neutrinos with m ~ keV the transfer function is suppressed on scales A < 900 kpc, the scale of suppression for 
TO ~ keV sterile neutrinos produced by scalar decay at a scale ~ 100 GeV is A < 488 kpc with a large enhancement of 
power at smaller scales. 

For sterile neutrinos produced by scalar decay we find the following behavior for the transfer function: at long 
wavelengths, 
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Tik)^l-Ci^-^-^^j +••• ; k^kMte,) (1.1) 

with C ~ 0(1), and at small scales where the corrections to the fluid description and the memory of gravitational 
clustering becomes important 

r(fc)~ 1.902 e-'=/'^'^=(*=«) ; k>kfs{teq) (1.2) 

valid for scales 65 kpc ^ A < 500 kpc where kfs{teq) is the free streaming wavevector at matter radiation equality. For 
TO ^ keV and decouphng temperature ~ 100 GeV we obtain kfs{teq) ~ 0.013/kpc. 

The smaller suppression scale may relieve the tension between the X-ray [3a| and Lyman-a forest [3y, [STj data and 
may provide the necessary enhancement of power at small scales to smooth out the inner profile of (dSphs). 

II. THE MODEL 

We study the model presented in references [23, [2^, [30, [SJ] as an extension of the minimal standard model with 
only one sterile neutrino, however, including more species is straightforward. The Lagrangian density is given by 
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C = CsM + -d,,xd^x " — X' + i^^^ - ^X^'^ " Y ^'^ " 2/ai^^^a i' - V{H^H- x) + h.c. (2.1) 

where Csm is the standard model Lagrangian, La]a = 1,2,3 are the standard model SU{2) lepton doublets, i/ is a 
singlet sterile neutrino with a (Majorana) mass m, a real scalar x with Yukawa coupling Y to the sterile neutrino 
which in turn is Yukawa coupled to the active neutrinos via the Higgs doublet -ff, thereby building a see-saw mass 
matrix in terms of the vacuum expectation of this doublet [27l. l29l. |34| . 

As discussed in detail in ref. |33j abundance and phase space density constraints from (dSphs) indicate that the mass 
of suitable (WDM) candidates must be in the keV range, which leads to considering_^the vacuum expectation value 
and mass M of the singlet scalar x in the range ~ 100 GeV as discussed in references [2^ |30, ISJ] • If the sterile neutrino 
mass m ~ keV results from the vacuum expectation value (x) ~ 100 GeV then the Yukawa coupling Y ^ 10"^. 

The results of the study here demonstrate that this range of parameters yields a consistent description of sterile 
neutrinos as a suitable (DM) candidate in this model. 

A. Decoupling out of equilibrium 



Non-resonant active-sterile mixing leads to sterile neutrino production via the Dodelson- Widrow mechanism 24 
with a decoupling temperature near the QCD scale t2^,i26|| ~ 150 MeV. In this scenario the distribution function at 
decoupling is of the form|24| 

U^{Pf-T) = -J^ ■ y = ^ ■ 0</3<l (2.2) 



where Pf is the physical momentum. For sterile neutrinos produced by non-resonant mixing with active neutrinos 24 
/3 ex 6*^ where 9m ^ 10~^[2a.l35J is the mixing angle. A fermionic relic decoupled in local thermodynamic equilibrium 
(LTE) while relativistic corresponds to /3 = 1. This general type of distribution function with a suppression factor /? 
has been used in the Lyman-a forest analysis ^STj. 

In this article we will neglect this production mechanism and focus on the production via the decay of the gauge 
singlet scalar field x which, as discussed in references [29l [30, |3J] lead to "colder" relics. However, we will compare 
the clustering properties of the distribution obtained via this mechanism and that of relics that decoupled with 
the generalized distributions (|2.2|) . postponing to another study the complete kinetic description that accounts for 
both pro cesses. Since these production mechanisms are effective at widely different scales {^ 100 GeV for x ~^ vv- 
decay[23,|30,l3J|, vs. ~ 150 MeV for Dodelson- Widrowf24i l2q) we expect that possible corrections from mixing will 
be subleading and certainly so for the small y region of interest. A more complete study is forthcoming. 

We consider the case in which the Yukawa coupling y ^ 1 and M :^ m and assume that the scalar field x is 
strongly coupled to the plasma and is in (LTE) with a Bose-Einstein distribution function^ 
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where fc is a comoving wavevector and 



T{t) = ^ (2.4) 

a{t) 



where Tq would be the temperature of the plasma today. During the radiation dominated era the Hubble expansion 
rate is given bv|40| 

Hit)^lMgHt)^ (2.5) 

where g{t) is the effective number of relativistic degrees of freedom. 



The case where the scalar is out of equilibrium has been considered in[29l |3l 



The details leading to the quantum kinetic equation for the production of sterile neutrinos via the decay of the 
scalar field x ^^^^ given in the appendix and the final result is given by eqn. (IA8I) . 

For Y^ ^ 1 we expect that neutrinos will decouple early and their distribution function will freeze-out with 
Up, rip <C 1. This expectation will be confirmed below self-consistently from the solution of the kinetic equation. 
Neglecting the neutrino population buildup in the kinetic equation (jASp . namely setting Up = Uq = 0, neglecting 
terms of order vi? /M"^ <C 1, taking the scalar field to be in LTE and replacing the momenta in (|A8p by their physical 
values, we find from 
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where Pf{t) = p/a{t) is the physical momentum, p is the comoving momentum and 
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where q± (t) are given by eqn. (|A9p in the appendix in terms of the corresponding physical momenta. These values 
are determined by the kinematic thresholds for scalar decay. 

We anticipate self-consistently, that for Y <^ 1 neutrinos decouple at temperatures Td^ m, namely when they are 
still relativistic, therefore we can safely neglect terms of order m? /T'^{t) < rri^/T^ ^ 1. Under this assumption (to 
be confirmed self-consistently below), using Pf{t)/T{t) = p/Tq and neglecting terms of order ni^ /M"^ ^ 1, the kinetic 
equation above simplifies considerably. It proves convenient to use the dimensionless variables 
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leading to the following form of the quantum kinetic equation. 
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where we have introduced 



A(r) = 



Y-2 



M 
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(2.10) 
(2.11) 

(2.12) 

(2.13) 



and the effective number of relativistic degrees of freedom depends on time through the temperature. 

Under the assumption M ^ to, for example taking M ~ 100 GeV, m ^ keV, the exponential in the numerator 
inside the logarithm can be neglected for all y ^ m? /A-P ~ 10^^^. Although we are interested in the small momentum 
region of the distribution function (small y) , the phase space suppression for small momentum in the integrals of the 
distribution function entails that we can safely neglect the contributions of such small values of y. This argument 
will be confirmed explicitly below. Therefore we safely neglect the numerator inside the logarithm in (j2.12p . In 
order integrate the rate equation (|2.12p we must furnish g{t). In the Standard Model this function is approximately 
constant in large intervals and features sharp variations in the regions of temperature when relativistic degrees of 
freedom either decay, annihilate or become non- relativistic [40'|. 

We will assume that g(t) is approximately constant in the region of (large) temperature before and during decou- 
pling, replacing the value of g by its average g over the range in which the rate (|2.12p is appreciable. The numerical 
analysis presented below justifies this assumption for a wide range of y. Therefore we replace 
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With these appr oximat ions and assumptions, eqn. (|2.12p is equivalent to the quantum kinetic equation obtained 
references [271 [29l |42| and can be integrated exactly. We obtain, 
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and r[a, 6] is the incomplete Gamma function. 

The rate (|2.12p vanishes as r — > 0, reaches a maximum and falls-off exponentially as t jAy — > oo. The maximum 
rate is larger for smaller values of y, this feature translates into an enhancement of the distribution function at small 
momenta, which will be at the heart of the important aspects of clustering studied in section pvp . 

The asymptotic behavior of the distribution function (|2.15p for r/Ay 3> 1 is given by 
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Figure ([T]) displays the rate (|2.12p (left panel) and the distribution function (|2.15p (right panel) as a function of r 
for several values of y. For y ^ 1 the distribution function is largest, reaching its asymptotic value for t < 1, whereas 
for large values y 3> 1 the distribution function is strongly suppressed and reaches its asymptotic form much later 
(see figUl). 

Hence, for small y, the region of distribution function most relevant for small scale structure formation [47|, decou- 
pling occurs fairly fast, at a decoupling temperature ^ M ^ m thus justifying the assumption of a constant g{t) and 
Tn/Td <C 1 with Td = Tq. 

As discussed above and in ref. [47| , the dark matter transfer function depends more sensitively on the small momen- 
tum (small y) region of the distribution function, the above analysis shows that for y < 1 the distribution function 
freezes out on a "time" scale t ~ 1, namely a temperature scale T(i/) ~ M, where t/ is the "freeze-out" time. 
Hence if the effective number of relativistic degrees of freedom g{t) varies smoothly within the temperature range in 
which the population freezes-out ^ M the approximation (|2.14p is justified and reliable. Therefore we take the value 
'g as the effective number of relativistic degrees of freedom at decoupling, for example for a freeze-out temperature 
Titf) - M - 100 GeV in the standard model g - 100 40]. 
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FIG. 1: Left panel: dn{y; r)/drA for y — 0.2, 0.5, 1. Right panel n{y; r)/A for the same values of y 



The distribution function at freeze-out is obtained by taking the r/Ay 
introduce the distribution function at decoupling 
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oo limit in eqn. (j2.15p . therefore we 
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^ Surprisingly the function gs (z) also determines the equation of state of the ideal non-relativistic Bose gas. 



This is the unperturbed distribution function that enters in the Boltzmann-Vlasov equation that determines the 
evolution of density and gravitational perturbations, and is displayed in fig. Q. It is remarkable that for y ^ 1 
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in striking contrast to the thermal Fermi-Dirac distribution function and to the one obtained from the (DW) mechanism 
proposed in refs.[2J, [25|, and closer to that of a (non-condensed) bosonic massless particle for which the distribution 
function at small momentum is ^ 1/y. 




FIG. 2; The distribution function fo{y) — n{y;oo). 



The divergence of the distribution function /o(j/) for y — > must be interpreted with care because we have neglected 
the build-up of the neutrino population in the quantum kinetic equation, furthermore we have also neglected terms 
~ m/T{tf) ^ m/M in the derivation. Neglecting the neutrino population in the kinetic equation requires that 
/o(y) <C 1 for all values of y, and neglecting the ratio m/T{tf) in the frequencies requires that y ^ m/T{tf) ~ m/M . 
These constraints imply a range of couplings for which the approximations leading to the final result of the distribution 
function are reliable. For example, taking m ~ keV ; M ~ 100 GeV neglecting the term m/T{tf) ~ m/M requires 
that y 3> 10^*. To obtain an estimate of the range in which the population build-up can be neglected, it is convenient 
to write 



A ~ 3 (y X 10^)M 1?^ 



100 GcV 

M 



(2.20) 



As discussed above (see also ref.[2^, |3J|) taking the expectation value of the scalar (x) ^ M ^ 100 GeV and assuming 
that the neutrinos obtain their mass via the Yukawa coupling with m ^ keV, this implies Y ^ 10^^ leading to the 
illustrative estimate A ~ 0.03. The condition for negligible neutrino population becomes 
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leading to y ::^ 10"''^ for which we can safely ignore m/T{tf) and the build-up of the neutrino population. This 
constraint may be implemented by introducing an infrared cutoff j/c — 10~ in the y— integrals of the distribution 
function, however, our analysis below shows that this is a mild constraint because of the phase space suppression in 
these integrals, hence the lower limit can be safely taken to y = 0. Values of 5 > 100 and M > 100 GeV yield smaller 
values of A and larger region of reliability for y. 

III. CONSTRAINTS FROM (DM) ABUNDANCE AND (DSPHS) PHASE SPACE DENSITY. 



The number density today of this DM candidate is 
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where 
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AC(5) ; C(5) = 1.037- •• 



Since this DM particle is non-relativistic today, its energy density is given by 

Pom = m no . 
Entropy conservation [40| entails that 
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where T„, 



2.348 X 10 ^ eV is the temperature of the cosmic microwave background today, and gj^ = g is the 



effective number of relativistic degrees of freedom at decoupling. 

The condition that this DM candidate contributes a fraction < i'dm < 1 to the (DM) density yields the following 
upper bound on the mass[33| 



m< 2.695 eV 
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where g is the number of internal degrees of freedom of the particle. For a decoupling temperature '^ 100 GeV for 
whichf403 g ^ 100 and A ~ 0.01 — 0.1 this upper bound is in the keV range. 

For comparison, for a fermion relic that decoupled with the general distribution function (|2.2p the corresponding 
upper bound becomes 



m< 3.593f|i')eV 
\(3gJ 



(3.6) 



For sterile neutrinos produced via the (DW) mechanism, (3 is adjusted to satisfy this bound with a given mass and 
gd ~ 30 (although g^ varies rapidly near 150 MeV because of the QCD phase transition or crossover) and for a thermal 
relic (/3 — 1) with m ^ keV it follows that gd > 500, namely such thermal fermion candidate must decouple at a 
temperature much higher than the electroweak scale. 

As discussed in ref. [33i] a generalized Tremaine-Gunn[43| bound that yields a lower bound on the mass is obtained 
from the coarse grained primordial phase space density 
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where n{t) is the number density of non-relativistic particles and (Pf) is the average of the squared physical momentum 
in the decoupled distribution function. This quantity is a Liouville invariant in absence of gravitational perturbations 
when the particle has become non-relativistic [33J. The observable quantity is p/(T^ where p is the mass density and a 
the one dimensional velocity dispersion, therefore we define the primordial phase space densitv|33j 
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The distribution function (|2.18p yields 
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whereas for fermions that decoupled while relativistic with the generalized distribution (|2.2p V ^ q 0x 1.963 x 10~'^[33|. 
Since the phase space density only diminishes during the merger process (violent relaxation) [43, [Sll] a lower bound 
on the mass follows 1331. 



m > 



[62.36 eV] 



W 



_4 p (km/s)^ 
cr3 Mo/kpc^ 



(3.11) 



The compilation of photometric data from fdSphs)|14| yields [■ ■ ■]^ ~ 1 — 2, taking the middle of this range as an 
estimate, the mass of the (DM) particle is bound in the region between the lower bound (|3.1ip and the upper bound 
p.Sp . namely 



H^<^<1.33 4 eV 

Taking as representative values g ^ 100; A ^ 0.05; g = 2 yields 

560 eV < m < 1330 eV 



(3.12) 



(3.13) 



constraing the mass in a fairly narrow window within the keV range. A similar conclusion (based on a similar analysis) 
was obtained in ref. [30j. 

Thus we see the consistency between the model with scales (x) ~ M ~ 100 GeV; Y ^^ 10~*; m ~ keV and the 
constraints from abundance and phase space density of (dSphs). 



IV. NON-RELATIVISTIC BOLTZMANN EQUATION: 

SPECTRUM 



TRANSFER FUNCTION AND POWER 



The transfer function is obtained from the solution of the linearized Boltzmann equation for density and gravitational 
perturbations. When the particle has become non-relativistic and for wavelengths that are well inside the Hubble 
radius, the non-relativistic Boltzmann- Vlasov equation describes the evolution of these perturbations. This equation 
was used in pioneering work on non-relativistic dark matter [5J|, for neutrinos [23, l33|, dark matter perturbations 
accreted by cosmic strings |55i [56| and more recently to study clustering of thermal relic neutrinos [57]. In all of these 
previous treatments a numerical analysis was offered but always with a thermal distribution function that is truncated 
to facilitate the numerical integration. 

Instead, here we follow the analysis of ref. [47| and implement very accurate analytic approximations for the transfer 
function that yield a deeper understanding of the connection between the decoupled distribution function and the 
transfer function, and provide a simple framework to obtain a reliable assessment of the power spectrum for arbitrary 
range of parameters (mass, coupling, etc) and distribution functions. 

To linear order in perturbations the distribution function of the decoupled particle and the Newtonian gravitational 
potential are [4^ [52, [53| 



f{p;x;t) = fo{p) + Fi{p;x-t) 

ip{x,t) = LpQ{x,t) + Lpi{x,t), 



(4.1) 
(4.2) 



where foip) is the unperturbed distribution function of the decoupled particle, given by (|2.18p and for comparison we 
will also study the generalized distribution (|2.2p with y = p/Tq, (po{x, t) is the background gravitational potential that 
determines the homogeneous and iso trop ic Friedmann-Robertson- Walker metric and p, x are comoving variables. The 
reader is referred to refs.[33, [4al53,|5J, [5^, [5g| for details on the linearization of the coUisionless Boltzmann- Vlasov 
equation. 
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In conformal time r and in terms of comoving variables p, x it is given by[33, 0, [53| 

IdFi p 



„ . 2 ^^^1 - mVs^i ■ Vpfo = (4.3) 

a OT ma'' 

along with Poisson's equation 

V|^i / — ^Fi(x,r) . (4.4) 

a J (27r)'^ 

It is convenient |54l l5a| to introduce a new "time" variable s related to conformal time r by 

dr 
ds=—, (4.5) 

a 

and to take spatial Fourier transforms of ipi {x, r) and Fi {x, r) obtaining 

^^'^I'P'^'^ + '±lF,{k,p;s)-ik-^j,fo{p)a^{s)Mk,s) = 0, (4.6) 

OS m 

.J* , 47rGm f d^p „ ,r ^ n /. -T^ 

The solution of equation (|4.6p is 

Fi(fc,p;s)=Fi(fc,p;sOe-'^(^^^')+imfc-Vp-/o(p) / ds' e^'^^^'-^-'a^^^/^^^^^^ ^.^ ^ (4 8) 



where 



The first term on the right hand side is the free-streaming solution in absence of gravitational perturbations. The 
analysis in ref.^3 shows that {p/m){s — si) is the free streaming distance that the particle travels with comoving 
velocity pim from Si until s. Multiplying both sides of eqn. (|4.8p by —A'i:Gml[k'^a{s)], integrating in p, and using 
the relation (|4.7p . we obtain 



Ba{s) J {2Tif'^ "pj^'yr, j^-- - ^ ^^ ,^.y^,^ , fc2a(s)7 (27r)3^ 

(4.9) 
The inhomogeneity on the right hand side of this equation is determined by the first term in (|4.8p and describes the 
free streaming solution of the Boltzmann-equation in absence of gravitational perturbations. 

During matter domination and choosing the initial time at matter-radiation equality Si — Seq, it is found that 47l| 

2u 

S-Seq^ , (4.10) 

HoMO-eq 



where 



^OM = — 5— POA/ = HqQm , (4-11) 



with Hq — 100 /i Km sec ^ Mpc is the Hubble parameter today and poM is the matter density today. 
Normalizing the scale factor to unity today, namely a(0) — 1, the variable u is given by[47| 



a / 


= 1- 


" 1 + Z ' 
1 + Zeq 


2 


nby 


^^„A 


^eq 



u=l-i-^\ =1- -;— ; 0<u<l-aiq, (4.12) 

and the scale factor in terms of u is given by 



' ' (l-u)2 ■ 
The initial value of the scale factor at matter-radiation equality is a^q — 1/(1 + Zeg) with z^q ~ 3050. 



(4.13) 
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It is convenient to introduce the normalized unperturbed distribution function 



My) = ,00 ^ff] ,, , (4.14) 

Jo y fo[y)dy 



which for the non-equihbrium distribution ()2.18p is given by 



4 95 (y) 
/o(;/)= '\' . (4.15) 

3V7rC(5) y2 



Following ref.([47|) we introduce the density perturbation normalized at the initial time 



J j2^Mk, P;Seq) 



and the gravitational perturbation normalized at the initial time 



with the relation (see eqn. (|4.7p ) 



*(fc,") = ^^ (4.17) 



$(fc; u) = -T^Sik; u) = (1 - uf6(k; u) . (4.18) 

a(u) 



The gravitational perturbation $ obeys Gilbert's equation[54| 

$(fc, u) - -(1 - u? f n\a(u - u')] P^^^ du' = (1 - ufllk, u] . (4.19) 

a Jq [1 - u'J* 



where we have introduced [41 



2fcl^l /100\' /keV 



AIlJ-^ ~ 1.18 k X [kpc] ( ^ 1 ( ^^^ 1 yiT^. (4.20) 



The non-local kernel is given by 



and 



n[z]= / dyyUy)sm[yz], (4.21) 

Jo 



ir P^'^P -^1 (^' P ; seg) Jo ( ^ ) 

/[fc, m] = ^^^ (4.22) 

/o p'^dpFi(k,p;Seq) 

where jo{x) — sin(a;)/a; and we have assumed that Fi is a function of \p\. The inhomogeneity I[k, u] is recognized as 
the free streaming solution normalized at the initial time, it obeys the initial conditions 

I[k, u = 0] = 1 : —I[k, u] = . (4.23) 

du M=o 

The density perturbation 5 obeys 

5[k,u)~- I main - u')] f^^' "^ du ^ I\k, u] . (4.24) 

a Jo [1 - "T 

In ref. 47] it is proven that 6 oc a(t) as w — > 1, just like fluid density perturbations in a matter dominated cosmology. 
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As analyzed in ref. [43|, the transfer function is obtained from the solution of Gilbert's equation (|4. 19^ ^47]. Normal- 
izing it so that T{k = 0) = 1, it is given bv[47t 



$(0,u) 3 
The final power spectrum Pf{k) is related to the initial one Pi{k) as[45| 

Pfik)^T^{k)P,{k). 



(4.25) 



(4.26) 



If perturbations do not grow or decay substantially during the prior, radiation dominated phase, Pi{k) is nearly the 
inflationary primordial power spectrumJ45|, which is given 



P^{k)^A\ 



ko 



(4.27) 



where the amplitude A and index n^ are determined during inflation and fcg is a pivot scale. The five years data 
release from WMAPfH] yields Us « 0.96. 

Initial conditions: The inhomogeneity I[k, u] depends on the initial condition Fi{k,p; Seq) which must be specified 
from the evolution of perturbations during the radiation dominated (RD) era up to t^q. During (RD) the Boltzmann 
equation (|4.3|) describes the evolution of DM perturbations when the particles are non-relativistic, for the cases under 
consideration T < m ^ IkeV, corresponding to a > lO^^aeq, and for gravitational perturbations with modes well 
inside the horizon. During this stage the evolution of the gravitational potential is determined by its coupling to the 
radiation fluid, it is given by ipi{k;t) = •i^pji(x)/x where $pis the amplitude of primordial perturbations, j'l is the 
spherical Bessel function, x = krj/^/i and 77 conformal time[45l|. ipi{k;t) is strongly suppressed for wavelengths well 
inside the horizon and can be considered a small perturbation to the evolution of (DM) perturbations [45]. Therefore 
the evolution of non-relativistic perturbations during (RD) is obtained from the solution (|4.8p by neglecting the 
second term that includes the gravitational potential, namely it is given solely by the free-streaming contribution. 
Consider an initial condition determined early in the (RD) era at s* . The free streaming solution of (|4.8p during (RD) 
(neglecting the gravitational potential) is 



Fi(fc,p; s) = F,ik,p; s*) e^'^^^^"^*) 



In 



a-(s) 



a(s') 



[Hq ^Im a-eq 



1/2 



(4.28) 



Therefore extrapolating this solution to matter-radiation equality, it follows that the initial condition at t^q, namely 

Fi{k,p;Seq) is given by (|4.28p evaluated at a{s) = Qeq- 

This analysis is akin to the evolution of (CDM) perturbations after decoupling from the radiation fluid studied in 
ref.[59|. An important difference however, is that sterile neutrinos cannot be described as part of a radiation fluid 
because they do not interact with radiation, leptons or baryons. 

We can now follow all the steps described above leading to eqns. (|4.19l4.24p and obtain the same equations but the 
inhomogeneity /[fc;u] replaced by 



I[k; 



I[k;u + uq] 
I[k;uo] 



Wo = 2 1^ 



^eq 



a{s* 



>0. 



where now 



I[k,u] = 



l^p^dpF,ik,p;s*)jo{^) 



Jo p'^dpFi{k,p;s*) 
taking Fi to be a function of \p\ . It is clear that defining the rescaled perturbations 

$(fc,u) = ^{k,u)I[k;uo] 
5{k,u) — 5{k,u)I[k]Uo] , 



(4.29) 



(4.30) 



(4.31) 
(4.32) 



these obey equations (|4.19l4.24p with the inhomogeneity I[k\u -f uq] where I[k\u\ is given by (|4.30p . As it will be 
shown below I[k; ug] < 1 for mq > 0. The rescaling of the gravitational and density perturbations by the factor /[fc; mq] 
reflects the suppression of perturbations by free streaming during the prior (RD) era. 
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It now remains to determine Fi{k,p; s*) to finally obtain I[k;u\, here we assume adiabatic initial conditions 
corresponding to a temperature perturbation, namely 



m,..>')=(T'J^){^) 



(4.33) 



These are the initial conditions proposed and studied in refs.[2l|, |39|. 

Although a more detailed analysis of the evolution during the radiation era and its impact on the small scale structure 
will be presented elsewhere, we can obtain an upper bound on the small scale properties of the transfer function by 
setting uq = {I[k; 0] = 1). Furthermore, since our objective is to compare the small scale properties of the transfer 
function for sterile neutrinos decoupled with the distributions (|2.18p and (|2.2|) . it is clear that the distribution (|2.18|) 
leads to a smaller suppression during (RD) because it favors the small momentum region, therefore yields a smaller 
free streaming velocity and a smaller free streaming length as compared to the distribution (|2.2p . 

Hence from the above discussion we conclude the following: i) setting uq — leads to an upper bound on T(k), ii) 
the distribution function (|2.2p leads to a larger free streaming suppression during (RD) as compared to the case of 
the distribution function (|2.18p . 

With the distribution function (|2.18p and initial condition (|4.33p we find that the normalized free streaming solution 
is given by 



I[k,u] 



1 



9\/2C(5) 



E 



1 



pn 



1 + ^ 



n + p 



a-" 



{auY 



n 
1 + 6^ 



p'^ = 71^ + iauf (4.34) 



Its asymptotic behavior for a u ^ 1 is 



I\z] 



au . 



(4.35) 



This power law fall off is strikingly different from that of weakly interacting massive particles (WIMPs) for which the 
fall-off by free streaming is exponential |47|. and that for fermions that decoupled with generalized distributions (|2.2p 
for which the normalized free streaming solution is the same as for a thermal relativistic relic (since the suppression 
factor (3 cancels out) and is also a power law but with a faster fall off oc z~'^, whereas for thermal bosons that 
decoupled in equilibrium the fall-off is ex z~^[47|. Thus the non-equilibrium distribution function fn(i/) p.l8p leads to 
a suppression of the free streaming solution with a power law intermediate between that of thermal ultrarelativistic 
fermions and bosons. Fig. ^ displays I[z];z^^^ I[z] vs. z = au. 




5/2 J. , 
Z I[z] 



FIG. 3: The free streaming solution I[z] eqn. (|4.34|l and z2l[z] vs. z — au. 



It is illustrative to compare the free-streaming solution (|4.34p to that of a fermionic species of the same mass that 
decoupled at the same temperature, therefore has the same a, and initial condition (|4.33p but with the generalized 
distribution (|2.2p . The free-streaming solution is independent of /3 and is the same as for a thermal relativistic relic [47| 



ILTE 



fc, u] — 



9C(3) 



y^ei 



sin[yau\ 
{ey + 1)2 an 



dy = 



3C(3) 



oo 

E 



(_1)("+1) 



+ Z^Y 



1-1 

3 



[n^ + z^ 



(4.36) 
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FIG. 4: The free streaming solution I[z] for the non-equilibrium distribution (|2.18p eqn. (|4.34p (sohd hne) and for the 
generahzed distribution function eqn. (|2.2|l (dashed hne) vs. z = au. 



It is clear that, at least for the initial conditions corresponding to temperature fluctuations (adiabatic) (|4.33p . the 
free-streaming solution in absence of gravitational perturbations has a slower fall-off in the case of the non-equilibrium 
distribution function when compared to the case of the generalized distribution (|2.2p . 

This remarkable difference also emerges in the non-local kernel Il[z] in Gilbert's equations for gravitational pertur- 
bations (|4T9l) or for density perturbations (|4.24|) . We find 



n[z] 



V2. 



V3C(5) ,§ 



[pnj 



2n- 



n + p 



\/ v? + z^ 



(4.37) 



For z ~ n[z] oc z and asymptotically for z — a{u — u') 3> 1 we find n[z] ex z ■^'^, in contrast with the case of thermal 
fermions for which the asymptotic behavior is Tl[z] ex z"^, whereas the asymptotic behavior for thermal bosons is 
found to be Il[z] oc ^"^[ig. Fig.® displays U[z] vs. z. 




FIG. 5: The kernel n[z] vs. z = a{u- u) 



As discussed in detail in ref. [47j the longer range of a kernel leads to an enhancement of the transfer function, and 
to more power at small scales. 

The free streaming wave vector: 

The comoving free-streaming wave vector is akin to the comoving Jeans wavevector in a fluid but with the speed 
of sound replaced by the velocity dispersion of the decoupled particle, it is given by 



kfs{t) = kfMVW) 



(4.38) 
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and its value today is|47j 



kMo) 

where 



2(y2) 



(4.39) 



(^'> = r}r/r , = (-) y' (4-40) 



is the three dimensional velocity dispersion of the non-relativistic particles today and we introduced 

/•OO 

F= / dyy'Uy). (4.41) 



Jo 

^2 



Using the relation (|3.4p and Jlji/Zi^ = 0.105 for non-baryonic DM[58[, leads to[4j 



,„, 0.563 / qd\3 / m \ ,_, . ^ ,, .^, 

fc/.(0) = ^(^) (i^)[kpc]-^ (4.42) 



The variable a defined in cqn. (|4.20p is related to kfs{ti) as 



6 \ k 

t kfs{teq) 



(4.43) 



For the non-equilibrium normalized distribution function (j4.15p we find 



whereas for the generalized distribution (|2.2p is the same as for thermal relativistic fermions, 

^Lt£ = 15 II ^12.939. (4.45) 

Therefore the enhancement of the non-equilibrium distribution function at small momenta yields a ~ 30% reduction 
in the squared velocity dispersion. Taking the number of relativistic degrees of freedom at decoupling gd = 9 (see 
eqn. (|2.14p and preceding discussion), we find 

.,.(O)^0aS3(I)'(j=L)h,cl-.0.71l(X)' (£L)|kp„,-. («6) 

leading to free streaming wavevector and wavelength at matter-radiation equality 

^/.(te,)^ 0.013(4)^ (j^)[kpc]- (4.47) 

. X 27r /100\^ /keV\ , , , , 

For comparison, the value of kf s{0) for a fermion that decoupled with the distribution function (j2.2p is the same 
as for a relativistic thermal fermion, given by [47i] 

4f) (0)^0.157 (I) ^(j|^)[kpc]-i. (4.49) 

The smaller velocity dispersion in the case of the non-equilibrium distribution function yields a ~ 30% increase in 
the free streaming wave vector and consequently a decrease in the free streaming length. But just as importantly in 
enhancing kfs{0) is decoupling at high temperature for a larger value oi gd- 



16 

For comparison, a sterile neutrino produced via the (DW) mechanism Td ^ 150MeV;.gd ^ 30[2J| yields a free 
streaming wavevector at matter-radiation equality /^(^ {teq) ~ 900 kpc. 

As discussed in ref. [431 ^^ analytic understanding of the transfer function is obtained by rewriting eqn. (|4.24p as 
a differential-integral equation. Taking two derivatives with respect to u and using the original Volterra eqn. (|4.24p . 
we obtain 

S{k, u) — J- ^-2■ + 37"5(fc, u) ~ du K[u — u] j- — ^,^ = I[k, u] + 3j I[k, u] , (4.50) 

[I ~ u) Jq [1 ~ u ) 

where dots refer to derivatives with respect to the variable u and the non-local kernel K[u — u'] is given by 



K[u~u']=&a / y{y^ ~ y^)h{v) sin[a y{u - u')] dy . (4.51) 

Jo 

and 7 is defined as 

372 ^a^P. (4.52) 

From the definitions (|4.2Qp and (|4.39p we recognize that the dimcnsionless ratio 

V2k 



kfs{teq) 



(4.53) 



It is convenient to write eqn. (|4.50p as 



6{k, u) - |ii^ + 3^2 ^(^^ ^) ^ g^^. ^] _ (4 54) 

(1 — uY 



The source term is 



S[5]u] = Sa[u] + Si[5]u], (4.55) 
where 

Sa[u] ^ i + il^I, (4.56) 

S^[5-u] = rdu'K[u-u']^^^. (4.57) 
Jo [^ ^ "^ \ 

Passing to cosmic time it is straightforward [47j to find that the left hand side of eqn. (|4.54p is precisely of the form 



of the Jeans equation for fluids but with the adiabatic speed of sound replaced by y (V^). However, as discussed 

in ref.[47[ the non-local kernel ^Ju — u'] includes higher order moments oi p^/m? than those typically kept in the 
hierarchy of moment equations [45|. 

The two terms in the source S"!!?;?/] have very different physical interpretations. The first term, Sq describes a 
"driving force" resulting from the free streaming of the initial perturbation, the second term 5*1 is a correction to 
the fluid description and can be interpreted as a non-local "pressure" term. As discussed in ref.[43| the second term 
is negligible in the long wavelength limit since K[u — u'] ex a^ for a — > 0, but becomes important at small scales. 
Furthermore the memory of this kernel is determined by the small- j; behavior of the distribution function /o {y) [47| : 
larger support at small values of y yields longer range memory kernels, which enhance the transfer function at small 
scales. This is a consequence of the fact that for kernels with longer range, memory of the initial stages of gravitational 
clustering persists throughout the evolution leading to a larger contribution to 5*1 [47| . 

The solution of (|4.54p can be written exactly in a formal iterative Fredholm series in which the main ingredients are 
the mode functions corresponding to the fundamental regular and irregular solutions of the homogeneous equation. 
Defining 

z = V37(l - u) = zo(l - w) ; zo = V37- (4.58) 

these are 

hi{z) = f — - 1 J cos(z) 4- - sin(z) , (4.59) 



h2{z) = ^ - 1 sin(z) - - cos(z) , (4.60) 



17 



with the foUowing asymptotic behavior as z — > (?i — > 1) 

3 



z 
15 



(4.61) 



We emphasize that the kernels n[z];_ftr[z] do not depend on the overall normalization of the distribution function. 
This is important because Dodelson-Widrow-type distribution functions (|2.2p (see also ref.[2^) are of the form 



fdwiv) 



a 



ey 



1 



where < /3 < 1 is a suppression factor. For these distribution functions the normalized counterpart 



fdwiv) 



fdiuiv) 



1 



iry'fdUy)dy ^m^y 



1 



(4.62) 



(4.63) 



is the same as that for a relativistic fermionic thermal relic. The suppression factor f3 only enters in the abundance 
and the primordial phase space density. Therefore for sterile neutrinos produced by the (DW) mechanism [24|, l2a| 
or for general distribution functions of the form (|2.2p . the kernels in Gilbert's equation (|4.19|) or in the equivalent 
equation (|4.50p and the free streaming solution in absence of gravitational perturbations I[k,u] are the same as for 
the case of a fermionic relativistic thermal relic. This results in that the transfer function and power spectrum for 
sterile neutrinos produced by non-resonant mixing are the same as that for relativistic thermal fermions (see below) . 
This observation is relevant in view of the stringent constraints from the analysis in ref.J37J|. The results in this 
reference are based on hydrodynamical simulations that assume a sterile neutrino distribution function of the form 
(|4.62p . The overall multiplicative normalization does not affect the non-local kernel K[z] and as a result, nor does it 
affect the transfer function and the power spectrurn (see below). This kernel features the asymptotic behavior oc 1/z^ 
for any distribution of the form (|2.2p as used in ref. tSTt] for the simulations, whereas for the distribution function (|2.18p 
K[z] ex 1/2:2. Namely, the kernel falls-off slower than in the case of the Dodelson-Widrow production mechanism 
leading to a longer range of memory of gravitational clustering. The analysis in ref.f47j indicates that this feature, 
slower fall-off and longer memory range translates in an enhanced transfer function at small scales. This behavior 
will be confirmed below. 



A. The transfer function and power spectrum 

From the Fredholm series solution for S the exact transfer function is given by[47| 



T{k) = 



10 



V3 



T JO 



h2{u') 



I[k,u'] Si[5-u'] 



[l-u'Y 



6 



du' . 



(4.64) 



where 6 in 5*1 is the Fredholm solution of the integral equation (|4.54p . As shown in detail in ref.[47[ a remarkably 
accurate approximation to the full transfer function is obtained from the first two terms in the Fredholm series 



T(fc)~TB(fc) + r(2)(fc) 

where the first. Born term is given by 

Tsik) 
and the second order correction is given by 

Ti2){k)^ 
where 6^^'{k;u) is given by 

5'^^\k,u)^I[k,u] + 



10 



1/37^ Jo 



h2iu') 



I[k,u'] 5i[<5(i);u'] 



[l-u')^ 



du' 



10 



, I[k,u \ , 
n2[u )— TTT^du 



V37370 ' '{i-wy 



h2iu')Si[S^-^^;u']du', 



3\/373 Jo 



\/37 Jo 



[hi{u)h2{u') - h2{u)hi{u')] — — '——rdu . 

(1 — u'Y 



(4.65) 



(4.66) 



(4.67) 



(4.68) 
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Since the free streaming solution I[k; u], the kernels and the mode functions are all functions of a (see eqn. (|4.20|) ). 
it follows that T{k) is a function of a. Therefore it proves convenient to relate the comoving wavelength A — 27r/fc to 
a in order to establish the scales that enter in T{k), namely 




(4.69) 

where gd = 9 for sterile neutrinos with the non-equilibrium distribution function (|2.18p . 

The contribution from Ti2\{k) is the first correction beyond the fluid approximation and includes the memory of 
the initial conditions and gravitational clustering. This correction is negligible in the long wavelength limit a — > 
but becomes important at short scales[47l]. 

Fig. ([6]) displays the Born term Tg{k), the second order corrected {Tsik) +T(2){k))^ and the exact T^{k) obtained 
from the numerical integration of Gilbert's equation (|4.19p with the distribution function (|2.18p . This figure shows 
that: (a) the second order correction becomes important at small scales a > 1 and (b) that the Born plus second 
order correction approximation to the transfer function is a remarkably accurate. The outstanding accuracy of the 
second order approximation was also pointed out in ref.[47| for thermal relics. 

|TB(k)+T2(k)|' 




3 4 5 

a 

FIG. 6: Comparison between the exact (red squares) solution, the Born approximation and the second order improvement. 

Fig. ([7]) compares r^(fc) for fermions that decoupled with the generalized distribution function (|2.2p . because the 
normalized distribution is the same as the case of thermal relativistic relics wc refer to this case as thermal, and with 
the non-equilibrium distribution function (|2.18D (non-equilibrium). The right panel displays ln(r^(fc)) for both cases 
to make explicit the enhancement of the transfer function for the non-equilibrium case at small scales a > 1. There 
are two different sources of this small scale enhancement: i) the initial condition has a slower fall-off with k in the 
non-equilibrium case and ii) the kernels 11 [z] and consequently K[z] also have a slower fall-off with k for large k. Both 
aspects are a consequence of the enhancement of the distribution function for small y. 

Distribution functions that favor the small momentum region yield memory kernels that fall off slower an enhance 
the transfer function and power spectrum at small scales [47]. This small scale (large a) enhancement is clearly 
exhibited in the comparison in fig. ([7]). 

Although the simple analytic approximation is easily to study numerically, it is illuminating to provide fitting 
formulae for T(fc) in different ranges of scales. For large scales k ^ kfs{teq) equation (|4.50p can be solved in 
perturbation theory in 7 (or a) and the contribution from the non-local kernel K can be neglected to leading order 
in a, since in the long wavelength limit K ex a^. In the long-wavelength limit k <C kfs{teq) we find that T(k) is 
approximated by 



with C ^ 0(1) and depends on the distribution function[47|. Of more interest is the small scale behavior for 
k > kfs{teq), because at small scales the contribution from the non-local kernel which is a correction to the fluid 
description that includes memory of gravitational clustering becomes important. Figure ^ shows that T(k) can be 
approximated by an exponential for a > 0.8. A numerical analysis yields 

T(fc)~ 1.902 e-^-/*'-^^^*-) ; fc > fc/,(te,) , (4.71) 
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FIG. 7: Left panel: comparison between T (fc) for thermal fermions and sterile neutrinos decoupled out of equilibrium, right 
panel: comparison between ln(T^(fc)) for both cases. 



in the range 0.8 < a < 6, where 



^/.M ^ 0.013 (X)^(^)[kpe]-, 



(4.72) 



The fit (|4.70p is better than 2% in this range. For g ^ 100 and ?7i ~ 1 keV the fit (|4.71|1 accurately describes the 
transfer function in the range 

65kpc< A< 5G0kpc, (4.73) 

approximately from the scale of clusters of galaxies to that of galaxies. Within this range 10^^ < T{k) < 0.7. 
The excellent fit (|4.7ip is different from the often quoted numerical fit by Bardeen etaL|60l|. 

Power spectrum: 

Since T(fc) is a function of the combination a given by eqn. (|4.20p it is convenient to write the power spectrum 



P(fc)=A(^| T\k) = Ba"' T^{k). 



(4.74) 



Fig. ([5]) displays P{k)/B vs. a for sterile neutrinos decoupled with (|2.2p (equilibrium) and (|2.18p (non-equilibrium) 
as a function of a. The figure reveals that at large scales a ^ 1 both feature the same power spectrum (which is 
the same as that for cold dark matter) but a substantial difference emerges at small scales. The non-equilibrium 
distribution function (|2.18p yields an enhanced transfer function, hence more power at small scales. This is a conse- 
quence of the fact that the non-equilibrium distribution function favors smaller values of the momenta (small values 
oiy), leading to smaller velocity dispersion hence effectively colder particles, smaller free streaming length, but more 
importantly a memory kernel of longer range. This feature results in that memory of the gravitational clustering 
"lingers" longer and the initial value of the gravitational potential influences the process of gravitational clustering 
during a longer period of time[47| leading to an enhancement of the transfer function and the power spectrum at 
small scales. 

Fig. d?]) shows that the suppression scale of T^{k) for rehcs that decoupled with (|2.2p is at fc ~ kfs{teq)- For 
sterile neutrinos produced via the Dodelson-Widrow mechanism at a temperature Td ~ 150 MeV, with gd ^ 30 and 
m = IkeV this scale corresponds to a comoving wavelength ^t^ (teq) ~ 0.9 Mpc, whereas for a keV sterile neutrino 
produced in the model under consideration that decoupled with (|2.18p with 5 ^ 100 this scale corresponds to a 
comoving wavelength Xfs{teq) ~ 0.49 kpc. At smaller scales, for a ^ 1 the difference in T'^{k) for thermal relics and 
the non-equilibrium distribution becomes more dramatic as shown in the right panel of fig. ([7]) . 

The small scale enhancement of T(fc) is a consequence of the small y behavior of the distribution function, which 
translates into a longer range kernel K\z\. Thus it is clear from these figures that the non-equilibrium distribution 
function, combined with the higher decoupling temperature, namely the colder behavior (smaller velocity dispersion) 
yield a substantial enhancement of power at small scales as compared to either thermal relics or to sterile neutrinos 
produced via non-resonant mixing with active neutrinos (namely a la Dodelson-Widrow). 
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FIG. 8: The spectrum P{k)/ B vs. a for the non-equihbrium distribution function (solid hne) compared to a thermal fermion 
relic (dashed line). 



For scales A ^ IMpc, namely a <^ 0.4 the transfer functions and power spectra of sterile neutrinos produced 
by scalar decay, the (DW) mechanism, relativistic thermal relics or (CDM) are essentially indistinguishable. The 
non-equilibrium case begins to feature larger power than (DW) and relativistic thermal relics for a > 0.5, namely for 
scales A < 0.8 Mpc and becomes substantially larger than either of these cases for small scales A < 490 kpc. 



V. CONCLUSIONS AND FURTHER QUESTIONS 



In this article we have implemented a program that begins with the microphysics of production and decoupling 
of a dark matter particle candidate, constrains the mass and couplings from the observed DM abundance and phase 
space density of DM dominated satellite galaxies (dSphs) and obtains the DM transfer function and power spectrum 
by solving the Boltzmann-Vlasov equation for density and gravitational perturbations during matter domination. 

The model studied is a phenomenologically appealing extension of the minimal standard model proposed in 
references [23, [23, [30, [SJ] in which sterile neutrinos are produced by the decay of a gauge singlet scalar. With 
the scale of the expectation value and mass of this scalar ^ 100 GeV a consistent description of ~ keV sterile neu- 
trinos decoupled strongly out of equilibrium at a decoupling temperature ~ 100 GeV emerges and satisfies the DM 
abundance and phase space constraints from (dSphs). 

The distribution function after decoupling for sterile neutrinos produced by gauge singlet decay features a strong 
enhancement at small comoving momentum ex l/^/p in contrast to sterile neutrinos produced via non-resonant mixing 
with active neutrinos (a la Dodelson-Widrow) for which the distribution function is that of ultrarelativistic thermal 
fermions multiplied by a suppression factor. Such distribution function was used in the hydrodynamical simulations 
performed in ref. [37| to analyze the Lyman-a forest data. 

We have implemented an accurate analytic approximation to the solution of the Boltzmann-Vlasov equation and 
the transfer function introduced in ref. [47| and obtained the power spectrum. This approximation allows to identify 
which features of the distribution function determine the small scale behavior of the transfer function. Distribution 
functions that favor small (comoving) momentum lead to longer range memory of gravitational clustering and slower 
fall-off of the free streaming solution. Both features lead to small scale enhancement of the transfer function and 
power spectrum. 

We compare the transfer function and power spectrum in the case of sterile neutrinos produced by gauge singlet 
decay and by non-resonant mixing with active neutrinos, and find that the former is substantially enhanced over latter 
at small scales. The transfer function and power spectrum of sterile neutrinos produced via non-resonant mixing is 
the same as that for fermionic thermal relics. The suppression factor in the distribution (sec cqn. 12. 2p modifies the 
abundance an primordial phase space densities but not the transfer function or power spectrum. 

While at large scales A ^ 1 Mpc the transfer function in both cases are nearly indistinguishable and the same as the 
(CDM) case, the power spectrum for m ~ keV (DW) sterile neutrinos produced by non-resonant mixing is suppressed 
below a scale A < 900 kpc whereas the transfer function for sterile neutrinos produced via scalar decay is suppressed 
below a scale A < 488 kpc and substantially enhanced at smaller scales when compared to the (DW) case. 
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We find the simple fits to T{k) in the hmits of large and small scales. For large scales we find 



2 



T(fc)^l-C(^^-^j +••• ; fc«fc;.(ie,), (5.1) 

with C '^ 0(1). In this long wavelength limit the fluid description is valid and the contribution from the memory 
kernel is subleading. 

At small scales the corrections to the fluid description in terms of the non-local kernel that includes memory of 
gravitational clustering becomes important, in the small scale regime k > kfs{teq) a simple and accurate numerical 
fit yields: 

r(/c)~ 1.902 e"''/'=^= (*='), (5.2) 

where kfs{teq) is the free streaming wave vector at matter-radiation equality. For a sterile neutrino with m ~ 1 keV 
decoupling at Td ^ 100 GeV we find 

fc/,(ieg)^0.013/kpc. (5.3) 

This fit is remarkably accurate in the wide range of scales 60 kpc ^ A < 500 kpc and is different from the often quoted 
result of ref . [60| ■ 

We have given arguments that show that the results presented above are an upper hound to the small scale properties 
of T{k), since the evolution of WDM perturbations during (RD) leads to further suppression of T(k) with a larger 
suppression for the case of sterile neutrinos with distribution functions of the form (|2.2p as compared to those with 
(|2.18[) . This is because the distribution function (|2.18[) favors the small momentum region leading to shorter free 
streaming lengths and larger free streaming wavevectors, allowing more power at small scales. A more detailed 
analysis of the initial conditions obtained by including the evolution during the (RD) will be reported elsewhere. 

The substantial difference between the suppression scales in the transfer function and power spectrum at small 
scales between sterile neutrinos produced by gauge singlet decay and those produced by the (DW) mechanism suggest 
that sterile neutrinos produced by scalar decay may relieve the tension between the constraints from X-ray [3^ and 
Lyman-a forest data[36ll37|. 

In order to assess whether sterile neutrinos produced by scalar decay may explain the cored profiles of (dSphs), a 
full N-body simulation with the power spectrum obtained above must be carried out. 

Whereas in this article we have focused on the production mechanism from gauge singlet decay, the model includes 
sterile-active mixing via a see-saw (Majorana) mass matrix. Therefore there is also a complementary mechanism of 
sterile neutrino production via active-sterile mixing akin to the (DW) mechanism, which is effective at a much lower 
temperature ^ 150 MeV. The wide separation of decoupling scales, ~ 100 GeV for scalar decay, vs. ~ 150 MeV for 
(DW) suggests that once the distribution function has been established after decoupling at the higher temperature, 
the non-equilibrium effects at the lower scale may not modify the small momentum region significantly. We conjecture 
this to be the case because the wide separation of scales suggests that the distribution function obtained from scalar 
decay, may be taken as an initial condition for the kinetics of production via active-sterile mixing, in which case for 
small mixing angle and neglecting again the build up of the population, the final distribution function would be the 
sum of ((TT5)) and ^^. 

We have not included this possibility in this article, postponing the more detailed kinetic description, a study of 
the consequences and different initial conditions to a forthcoming article. 

In obtaining the transfer function and power spectrum by solving the Boltzmann-Vlasov equation for (DM) and 
gravitational perturbations during matter domination, we have neglected the contribution from baryons and photons. 
Although these will only affect the results for T{k) and P{k) at the few percent level, as discussed in the introduction, a 
precise assessment of P{k) to few percent accuracy requires solving the full set of Boltzmann-equations by modifying 
publicly available codes. The study presented in this article provides a reliable preliminary assessment of (DM) 
candidates, allows a systematic comparison and highlights the important small-scale aspects, perhaps eventually 
justifying the more numerically demanding task of solving the full set of coupled Boltzmann equations for photons, 
baryons, gravitational and (DM) perturbations. 
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APPENDIX A: THE QUANTUM KINETIC EQUATION 

With the Yukawa mteraction £/ — Y x^v and x ^ scalar field with mass M and v either a Dirac or Majorana 
fermion field of mass m the quantum kinetic equation is obtained just as in Minkowski space time by obtaining the 
total transition probabilities per unit time of the decay and inverse decay processes. 

The quantum kinetic equation for the neutrino population n(p; t) (and similarly for the antineutrino) is of the usual 
form 



dn{p; t) 
Jt 



Gain — Loss 



(Al) 



where the gain and loss terms are obtained from the corresponding transition probabilities jA^^^ip. The processes are 
depicted in fig. ([9]). 





k 

X 



FIG. 9: The "gain" and "loss" contributions to the quantum kinetic equation from the decay x ^^ '^'^ ^^i^d inverse decay 
77i^ — » X) to lowest order in the Yukawa coupling. 

The gain term is obtained from the decay reaction % ^ I7jy depicted in the first term in fig. ^ corresponding to 
an initial state with Nk quanta of the scalar field x and fip^s^riq^si quanta of neutrinos and antineutrinos respectively 
with Nk — 1, f^p.s + 1, ng,5' + 1 quanta in the final state respectively. The corresponding Fock states are 



\i) ^ \Nk,np^s,nq^s') ; |/) = |iVfc - 1, np,^ + 1, n^,,,' + 1) 



(A2) 



The loss term is obtained from the inverse decay reaction V v —f x corresponding to an initial state with Nk quanta 
of the scalar field x ^^^id nps,nq^s' quanta of neutrinos and antineutrinos respectively with N^. — 1, Ups — l,7ig,s' — 1 
quanta in the final state respectively. The corresponding Fock states are 

\i) = \Nk,np^s,ng^s') ; \f) = \Nk + l,np,s ~ l,nq.s' - I) (A3) 

The calculation of the matrix elements is standard, the fields are quantized in a volume V in terms of Fock creation- 
annihilation operators, with the corresponding spinor solutions for the neutrino fields. To lowest order in the Yukawa 
coupling. 



M 



fi 



gain V V^ V 2r2fc 

similarly, for the loss term we obtain, 

Y h.p+n 



Y Sj: ^.^ _ 

''' '' ■ 'Nk^yi - Upy^l -Uq Ua{p,s)Va{q,s'){2'K)6{Vlk - UJp- UJq) 



M 



/'■ 



= —I- 

loss V V^ v 20fe 



^V^l + Nky/n^y^ Va{q, s')Ua{p, S) (27r)(5(r2fc ~ Up ~ Uq) 



(A4) 



(A5) 



where the spinors have been normalized to one and the frequencies 



rjfc = yjk^ + M 2 ; Up = Vp^ + 



7n' 



(A6) 



Summing the respective jAl/ip over k,q,s,s' and taking the infinite volume limit we obtain the total transition 
probability for the gain term 






gam 



^2 
47r 



= T ^— I a q — — ^— [LUpLLjq ~ p ■ q~ m J 



^ ^ftA-n 



p+q^p ^q 



Np+qi^ -"-p) (1-"?) 



(A7) 
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where T is the total reaction time. For the loss term we find the same expression but with the replacement Np^^ -^ 
1 + Np^g' , (1 — Up) ^ Up , (1 — fiq) ^ riq. Carrying out the angular integral using the energy-conserving delta 
function we obtain the final expression for the neutrino production rate 



dn{p;t) 1 



dt 



T 



[ E i-^/'i 



E 1-^ 



gain ^ ^ 

k,q,s,s' 



/»i 



loss. 






■?+ qdq 

UJr, 



Np+q i'^-'^p) (l-"g)-(l+%+<f) ^p Uq 

(A8) 



k^q^s.s' 

where q± are the roots of the equations 

{p ± (J±)^ + Af ^ = UJq^ + LJp 

As discussed in section ([11]) the relevant limit is M ^ m. In this limit M ^ m we find these roots to be 

Af2 



q± 



2m? 



(ujp ± p) 



The extension to the cosmological case replaces the momenta by the physical momenta 

P 



p-^Pfit) 



i{t) 



(A9) 



(AlO) 



(All) 



and for Major ana neutrinos n = n. 
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